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Abstract 

Several researchers proposed using non-Euclidean metrics on point 
sets in Euclidean space for clustering noisy data. Almost always, a 
distance function is desired that recognizes the closeness of the points 
in the same cluster, even if the Euclidean cluster diameter is large. 
Therefore, it is preferred to assign smaller costs to the paths that stay 
close to the input points. 

In this paper, we consider the most natural metric with this prop¬ 
erty, which we call the nearest neighbor metric. Given a point set P 
and a path 7, our metric charges each point of 7 with its distance to 
P. The total charge along 7 determines its nearest neighbor length, 
which is formally defined as the integral of the distance to the input 
points along the curve. We describe a (3 -I- £)-approximation algorithm 
and a (1 -f e)-approximation algorithm to compute the nearest neigh¬ 
bor metric. Both approximation algorithms work in near-linear time. 

The former uses shortest paths on a sparse graph using only the input 
points. The latter uses a sparse sample of the ambient space, to find 
good approximate geodesic paths. 
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1 Introduction 


Many problems lie at the interface of computational geometry, machine 
learning, and data analysis, inclnding-bnt not limited to: clnstering, man¬ 
ifold learning, geometric inference, and nonlinear dimensionality reduction. 
Although the input to these problems is often a Euclidean point cloud, a 
different distance measure may be more intrinsic to the data. In particnlar, 
we are interested in a distance that recognizes the closeness of two points in 
the same cluster, even if their Euclidean distance is large, and conversely, 
recognizes a large distance between points in different clusters, even if the 
Euclidean distance is small. For example, in Figure 1 , the distance between 
a and b must be larger than the distance between b and c. 
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b 
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Figure 1: The intrinsic density-based distance can differ from the Euclidean distance. 


There are at least two seemingly different approaches to define a non- 
Euclidean metric on a finite set of point on M'^. The first approach is to 
form a graph metric on the point set. An example of such a graph is the 
feth nearest neighbor graph where a point is only connected to another point 
if one is a kth nearest neighbor of the other. The edge weights may be a 
constant or the Euclidean distances. In this paper we consider the complete 
graph but the edge lengths are a power of their Euclidean lengths. We are 
particnlarly interested in the squared length, which we will refer to as the 
edge-sqnared metric. 

The second approach is to endow all of with a new metric. We start 
with a cost function c : —>■ M, which takes the point cloud into account. 
Then, the length of a path 7 is the integral of the cost function along the 
path. 


4(7) 



c(7(t)) 




dt. 


( 1 ) 


The distance between two points x, y G is then the length of the shortest 
path between them: 


dc{x,y) = inf 4(7). 
7 


( 2 ) 
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Note that the constant function, c{x) = 1 for all x G gives the Euclidean 
metric; whereas, other functions allow space to be stretched in various ways. 

In almost all applications mentioned above for cost-based metrics, in 
order to reinforce paths within clusters, one would like to assign smaller 
lengths to paths that stay close to the point cloud. Therefore, the simplest 
natural cost function on is the distance to the point cloud. More precisely, 
given a finite point set P the cost c{x) for x G is chosen to be N(x), the 
Euclidean distance from x to the nearest point in P. The nearest neighbor 
length (N-length) £n(7 ) of a curve is given by (1), where we set c(x) = N(x) 
for all points x G C. We refer to the corresponding metric given by (2) as 
the nearest neighbor metric or simply the N-distance. 

In this paper, we investigate approximation algorithms for N-distance 
computation. We describe a (3 -|- e)-approximation algorithm and a (1 -|- e)- 
approximation algorithm. The former comes from comparing the nearest 
neighbor metric with the edge-squared metric. The latter is a tighter approx¬ 
imation that samples the ambient space to find good approximate geodesics. 

1.1 Overview 

In Section 4, we describe a constant factor approximation algorithm ob¬ 
tained via an elegant reduction into the edge-squared metric introduced 
by [BRSll] and [VB03]. This metric is defined between pairs of points in 
P by considering the graph distance on a complete weighted graph, where 
the weight of each edge is the square of its Euclidean length. We show 
that the N-distance and edge-squared metric are equivalent up to a factor 
of three (after a scaling by a factor of four). As a result, because spanners 
for the edge-squared metric can be computed in nearly linear time [LSV06], 
we obtain a (3 -|- e)-approximation algorithm for computing N-distance. 

Theorem 1.1. Let P be a set of points in and let x,y € P. The nearest 
neighbor distance between x and y can be approximated within a (3 -|- e) 
factor in 0(n log n -|- time, for any 0 < e < 1. 

In Section 5, we describe a (1 -t- e)-approximation algorithm for the N- 
distance that works in time e~^^^^nlogn. Our algorithm computes a dis¬ 
cretization of the space for points that are sufficiently far from P. Never¬ 
theless, the sub-paths that are close to P are computed exactly. We can 
adapt our algorithm to work for any Lipschitz cost function that is bounded 
away from zero; thus, the algorithm can be applied to many of the scenarios 
describe in Appendix 2. 
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Theorem 1.2. For any finite set of points P C and any fixed number 
0 < e < 1, the shortest N-distance between any pair of points of the space 
can be (1 + e)-approximated in time 0{e~^^^'^n\ogn). 

2 Related Work 

Computing the distance between a pair of points with respect to a cost func¬ 
tion encompasses several significant problems that have been considered by 
different research communities for at least a few centuries. As early as 1696, 
Johann Bernoulli introduced the brachistochrone curve, the shortest path in 
the presence of gravity, as “an honest, challenging problem, whose possible 
solution will bestow fame and remain as a lasting monument” [Ber96]. With 
six solutions to his problem published just one year after it was posed, this 
event marked the birth of the held of calculus of variations. In this section, 
we review some work related to computing shortest paths in a weighted 
domain. 

Models for Geometries. The cost function c(-) in (1) can be deliberately 
picked so that the metric dc(-) coincides with well-known metrics. For ex¬ 
ample, c{x,y) = 1 gives the Euclidean metric in the plane and c{x,y) = Ify 
gives the Poincare metric in the half-plane model of hyperbolic (Lobachevsky) 
geometry [KatlO]. 

Motion Planning. Rowe and Ross [RR90] as well as Kime and Hes- 
panha [KH03] consider the problem of computing anisotropic shortest paths 
on a terrain. An anisotropic path cost takes into account the (possibly 
weighted) length of the path and the direction of travel. Note that this 
problem can be translated into the problem of computing a shortest path 
between two compact subspaces of under a certain cost function. 

Computational Geometry. Indeed, computational geometers are inter¬ 
ested in different versions of this problem. In the simplest case, c takes 
values from {l,oo}, i.e., the space is divided into free space and obstacles. 
This problem can be solved in polynomial time using visibility graph in two 
dimensions, and it can be e-approximated in three dimensions. For example, 
the computation of the Frechet distance can be posed in this way [AG95]. A 
slightly more complicated case occurs when we let c be a piecewise-constant 
function. Mitchell and Papadimitriou [MP91] formulated this problem in 
two dimensions and designed a linear-time algorithm to find the solution 
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within e-accuracy. They list the problem for more general cost functions 
as an open problem (See Section 10, problem number (3)). A series of 
works [ALMS98, AMSOO, RSOO, AMS05] has resulted in an e-approximation 
algorithm that computes the shortest paths in a weighted polyhedral surface 
in 0((n/\/e) log(n/e) log(l/e) time. 

Machine Learning. Sajama and Orlitsky [SO05] first applied density- 
based distance (DBDs) to semi-supervised learning. Assuming that the 
sample points are taken from an underlying distribution with density /, a 
density-based distance can be defined by setting c(x) = in (1) and (2).^ 

The goal here is to place points that can be connected through dense re¬ 
gions into the same cluster. Vincent and Bengio [VB03] and Bousquet et 
al. [BCH04] suggest estimating / using a KDE and then approximating the 
metric by discretizing the space in a similar fashion to Tsitsiklis [Tsi95]. 
However, they do not provide any analysis on the complexity of the dis¬ 
cretized space. Bijral et al. [BRSll] bypass estimating / by building a 
complete graph over a set of points {xi,X 2 ,.. • ,Xn} sampled with respect 
to /, in which the length of an edge {xi,Xj) is \ \xi — Xj\\p for fixed p and q, 
and computing pairwise shortest paths in this graph. Hwang et al. [HDI14] 
prove that for certain values of p and q, the latter metric and the density 
based metric are equivalent up to a linear factor for sufficiently large values 
of n. 

The nearest neighbor metric can be viewed as a special case of density- 
based distance when the underlying density is the nearest neighbor density 
estimator. 


3 Preliminaries 

In this section, we define some basic concepts that are used in the paper. 

3.1 Metrics 

In this paper, we consider three metrics. Each metric is defined by a length 
function on a set of paths between two points of the space. The distance 
between two points is the length of the shortest path between them. 

^The definition of DBDs in [SO05] is more general in that it allows for a choice of 
discount functions (we take the inverse raised to the power 1/d as this seems to be the 
most natural way of turning a measure of volumetric density to a measure of length). 
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Euclidean metric. This is the most natural metric defined by the Eu¬ 
clidean length. We use £( 7 ) to denote the Euclidean length of a curve 7 ; 
£( 7 ) can also be defined by setting c{x) = 1 for all x G in (1). We use 
d(x, y) to denote the distance between two points x, y G based on the 
Euclidean metric. 

Nearest neighbor metric. As mentioned above, the nearest neighbor 
length of a curve with respect to a set of points P, is defined by setting c(-) 
to be N(-) in (1). The nearest neighbor length of a curve 7 is denoted by 
^n( 7 )) and the distance between two points x,y with respect to the 

nearest neighbor metric is denoted by dN(x,y). 

Edge-squared metric. Finally, the edge-squared metric is defined as the 
shortest path metric on a complete graph on a point set P, where the length 
of each edge is its Euclidean length squared. The length of a path 7 in this 
graph is naturally the total length of its edges and it is denoted by 4 q( 7 )- 
The edge-squared distance between two points x, y G P is the length of the 
shortest path and is denoted by dsq(x,y). 

3.2 Voronoi Diagrams and Delaunay Triangulations 

Let P be a finite set of points, called sites, in for some n > 1. The 
Delaunay triangulation Del(P) is a decomposition of the plane into simplices 
such that for each simplex a G Del(P), the Delaunay empty circle property 
is satisfied; that is, there exists a circle C such that the vertices of a are on 
the boundary of C and int(C') n E is empty. The Voronoi diagram, denoted 
Vor(P), is the planar dual to Del(P). We define the in-ball of a Voronoi 
cell with site p to be the maximal ball centered at p that is contained in the 
cell. The inradius of a Voronoi cell is the radius of its in-ball. We refer the 
reader to [DBVKOSOO] for more details. 

4 Nearest Neighbor Distance Versus Edge-Squared 
Distance 

In this section, we show that the nearest neighbor distance of two points 
x,y £ P can be approximated within a factor of three by looking at their 
edge-squared distance. More precisely, dsq(x,y)/4 > dt^{x,y) > dsq(x,y)/12 
(see Lemma 4.1 and Lemma 4.4). 
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As a consequence, a constant factor approximation of the N-distance can 
be obtained via computing shortest paths on a weighted graph, in nearly- 
quadratic time. This approximation algorithm becomes more efficient, if the 
shortest paths are computed on a Euclidean spanner of the points, which is 
computable in nearly linear time [Hpll]. A result of Lukovszki et al. (Theo¬ 
rem 16(ii) of [LSV06]) confirms that a (l-)-e)-Euclidean spanner is a (l-|-e)^- 
spanner for the edge squared metric. Therefore, we obtain Theorem 1.1. 

Before, starting the technical part of this section, we remark that both 
the nearest neighbor and edge-squared metric can have 11 (log n) doubling 
dimension. An illustrative example is a star with dense points sets on its 
edges. 

4.1 The Upper Bound 

We show that the edge-squared distance between any pair of points x,y € P 
(with respect to the point set P) is always larger than four times the N- 
distance between x and y (with respect to P). To this end, we consider any 
shortest path with respect to the edge-squared measure and observe that its 
N-length is an upper bound on the N-distance between its endpoints. 

Lemma 4.1. Let P = {pi,P 2 , • ■ ■ ,Pn} be a set of points in and let 
Tn and dgq be the associated nearest neighbor and edge-squared distances, 
respectively. Then, for any distinct points x,y & P, we have that dN(x, y) < 

jdsq(x,7/). 

Proof. Consider the shortest path x = qi ^ q 2 ^ ^ Qk = y with respect 

to the edge-squared metric. Let 7 be the same path in parameterized by 
arc length that uses straight line segments between each pair qi, gi+i. By 
the definition of the edge-squared distance, we have 

k-l 

dsq(x,y) = '^d{qi,qi+if. 
i=l 

On the other hand, by the definition of N-distance we have 

dN(a;,y) < ^( 7 ) = [ N(s)fis. 

J-y 
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The following sequence of equalities follow by basic rules of integration. 


fc-i 


N(s)ds = 


N(s) ds 


j=l + l 

fc -1 / , 


< 


SU- 

£(/ 

i=i V<1.- 


fc-i 

E 

fc -1 

E 

2=1 

-d,q(x,2/). 


, li+H+l 


. <Ji+9i + l 
2 


N(s) ds + 


d(s,gi)(is 


2 + l 


N(s) ds 


' Qi+Qi+1 

2 ~^9i+l 


d(s,g'j+i)ds 


f-d(qi,qi+i)/2 f-d{qi,qi+i)/2 > 

/ tdt+ tdt 

'o Jo ! 


d{qi,qi+if 


Thus, the proof is complete. 


□ 


4.2 The Lower Bound 

Next, we show that the edge-squared distance between any pair of points 
from P cannot be larger than twelve times their N-distance. To this end, we 
break a shortest path of the N-distance into segments in a certain manner, 
and shadow the endpoints of each segment into their closest point of P to 
obtain a short edge-squared path. The following definition formalizes our 
method of discretizing paths. 

Definition 4.1. Let P = {pi,P 2 , • • • ,Pn} be a set of points in and let 
x,y ^ P. Let 7 : [0,1] —)■ be an {x, y)-path that is internally disjoint from 
P. A sequence Q < to < ti <■■■< tk < I is a proper breaking sequence 
of 7 if it has the following properties: 

1. The nearest neighbors of'y{to) and jitk) in P are x and y, respectively. 

2. For alll<i<k, we have £( 7 [fj_i, f*]) = ^{N{'y{ti-i)) + N{'y{ti))) 

The following lemma guarantees the existence of breaking sequences. 

Lemma 4.2. Let P = {pi,P2,-" :Pn} be a set of points in and let 
x,y G P. Let j be a path from x to y that is internally disjoint from P. 
There exists a proper breaking sequence of 



Proof. Pick to such that the closest neighbor to 7 (to) in P is x. Inductively, 
pick ti > ti-i so that property ( 2 ) of a proper breaking sequence holds until 
7 (ti)’s closest neighbor is y. 

We need to prove two properties: (I) a ti with property ( 2 ) always 
exists, ( 11 ) this process ends after a finite number of steps (i.e., ti falls in 
the Voronoi cell of y for some i). 

Property I. Let L_i be the last selected point in the process. We prove 
that a ti £ (L-i, 1 ) exists such that 

i{-f[ti-i,ti]) = ^(N(7 (L_i))+ N(7 (L))). 

Let the functions / : [L-i, 1] —)■ M and g : [U-i, 1] —)■ M be defined as 
follows: 

fit) = ^i7[ti-i,t]) 

ait) = ^(N( 7 (ti_i))+ N( 7 (t))). 

In particular, /(L-i) = 0 and g{ti-i) = N( 7 (ti_i)) > 0; so g{ti-i) > /(L-i). 
On the other hand, 

/(I) = £{-f{ti-i, 1 )) > N( 7 (ti_i)) > N( 7 (ti_i ))/2 = g{l). 

The first inequality is a result of the Lipschitz property and the fact that 
N(7(1)) = N( 2 /) = 0. 

Since both / and g are continuous functions, the intermediate value 
theorem implies the existence of a G such that /(L) = giU), 

which in turn implies property (I). 

Property II. For t £ [to, 1], let p{t) be the Euclidean distance from 7 (t) 
to P \ {y}. By definition, p{t) is positive everywhere. Since p is continuous 
and defined on a closed interval, by the extreme value theorem, it attains a 
minimum value Pmin > 0. In any inductive step, if the nearest neighbor of 
neither L_i nor L is y then i{'y[ti-i,ti]) > Pmin, by the second property of a 
breaking sequence. This implies the existence of a finite k such that 7 (tfe)’s 
nearest neighbor in P is y. □ 

Lemma 4.3. Let P = {pi,P 2 , ■ ■ ■ ,Pn} be a set of points in Furthermore, 
let 7 he any path in and x he an endpoint of 7 . If £( 7 ) = s, then 
^( 7 ) > sN(x) — s^/2. 
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Proof. Let be a unit speed reparameterization of 7 (i.e., | 7 u(i)| = 1 for 
all t). Suppose, without loss of generality, that x = 7 ( 0 ) = 7m(0). Then, by 
definition of the nearest neighbor metric, 

^n(7) == [ N{'yu{t))dt. 

Jo 

Then, the Lipschitz property of the N function implies: 


^n(7) > 


(N(7„(0)) - t)dt 


/ (N(x) — t)dt 

Jo 

sN(x) — s^/2. 


□ 

Given a path 7 that realizes the nearest neighbor distance between two 
points X and y, in the proof of the following lemma we show how to obtain 
another (x, y)-path with bounded edge-squared length. The proof heavily 
relies on the idea of breaking sequences. 

Lemma 4.4. Let P = {pi,P 2 ,--- ,Pn} be a set of points in and let 
dN and dgq be the associated nearest neighbor and edge-squared distances, 
respectively. Then, for any distinct points x,y € P, dj^{x,y) > ^dsq{x,y). 

Proof. Let 7 be a path from xtoy that realizes the nearest neighbor distance 
between x and y, i.e., ■^n( 7 ) = dN(x,y). 

Suppose, without loss of generality, that 7 intersects P only at {x,y}. 
Otherwise, we break 7 into pieces that are internally disjoint from P and 
prove the bound for each piece separately. 

Let ■ ■ • ,tk} be a proper breaking sequence of 7 . For 0 < i < A;, 

let Ui G P he the nearest neighbor of 7 (ti); in particular no = x and = y. 



Figure 2: Proof of Lemma 4.4 
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We show that > ^(d(nj_i, nj))^ for any I <i < k, which 

in turn implies 


%(7) > ^5^(d(ni_i,ni))2 > ^d^q{x,y). 

i=l 

Property 2 of a breaking sequence implies that = 2s, where 

■s = i(N( 7 (ii-i)) + N( 7 (ti))). We pick G so that i{'y[ti_i,mi]) = 

= s. Lemma 4.3 implies 

> (sN(7(ti_i)) - sV2) + (5N(7(L)) - sV2) 

= (3) 

On the other hand, by the triangle inequality, 

d(ni_i,rii) < N( 7 (L_i)) + 2s + N( 7 (ti)) = 6s. (4) 

The last equality follows by property I. Finally, Inequalities (3) and (4) 
imply 

^N(7i) > Y^(d(ni_i,ni))^ 

and the proof is complete. □ 

5 A (1+ £)-Approximation for the Nearest Neigh¬ 
bor Metric 

In this section, we describe a polynomial time approximation scheme to 
compute the N-distance between a pair of points from a finite set P C 
The running time of our algorithm is e~^^‘^^nlogn for n points in 
d-dimensional space. We start with Section 5.1, which describes an exact 
algorithm for the simple case in which P consists of just one site. Section 5.3 
describes how to obtain a piecewise linear path using infinitely many Steiner 
points. Section 5.4 combines ideas from 5.3 and 5.1 to cut down the required 
Steiner points to a hnite number. Finally, Section 5.5 describes how to 
generate the necessary Steiner points. 
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5.1 Nearest Neighbor Distance with One Site 

We describe a method for computing dN for the special case that P is a single 
point using complex analysis. This case will be important since distances 
will go to zero at an input point and thus we must be more careful at 
input points. Far away from input points, we will use a piecewise constant 
approximation for the nearest neighbor function but near input points we 
will us exact distances. More than likely this case has been solved by others 
since the solution is so elegant. We refer the interested reader to [Str] for 
more general methods to solve similar problems in the field of calculus of 
variations. 



Figure 3: To make the complex function one-to-one, one needs to extend the complex 
plane to the two-fold cover called the two-fold Riemann Surface. 

Suppose we want to compute dN(x,?/) where P = {(0,0)}. Writing 
{x,y) G C in polar coordinates as z = re*®, we define the quadratic transfor¬ 
mation f: C —)■ P by 

fiz) = zV 2 = (rV 2 )e* 2 ^ 

where TZ is the two-fold Riemann surface; see Figure 3. The important point 
here is that the image is a double covering of C. For example, the points 1 
and —1 are mapped to different copies of 1/2. Therefore, on the Riemann 
surface, the distance between 1 and —1 is one and the shortest path goes 
through the origin. More generally, given any two nonzero points p and q 
on the surface, the minimum angle between them (measured with respect 
to the origin) will be between 0 and 27r. Moreover, if this angle is > vr, then 
the shortest path between them will consist of the two straight lines \p, 0 ] 
and [0,g]. Otherwise, the line [p,q] will be a line on the surface and, thus, 
the geodesic from p to q. 

Let denote the distance on the Riemann surface. We next show that 
for a single point, the nearest neighbor geodesic is identical to the geodesic 
on the Riemann surface. 

Lemma 5.1. Let 7 : [0,1] —)• C 6 e a curve. Then, the image of 'y under f, 
denoted by f o j satisfies the following property: 

finif 01) = ^( t )- 
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Proof. Suppose 7 ; [0,1] —)• C is any piecewise differentiable curve, and let 
a := / o 7 . The N-length ^n( 7 ) of 7 is the finite sum of the N-length of all 
differentiable pieces of 7 . If the path 7 goes through the origin, we further 
break the path at the origin so that a is also differentiable. Thus, it suffices 
to consider (a, b) C [0,1] so that 7 [a, b] is a differentiable piece of 7 . Then, 
we have 

I • I is modulus. 

Modulus commutes with product. 
Chain rule. 


J a 

= f h{tW{t)\dt 

J a 

= j \a{t)\ dt 

J a 


□ 

Corollary 1 (Reduction to Euclidean Distances on a Riemann Surface). 
Given three points x, y, and p in such that p = NN(x) = NN(?/), the 
nearest neighbor geodesic G from x to y satisfies the following properties: 

1. G is in the plane determined by x,y^p. 

2. (a) If the angle formed by x,p,y is 7 r /2 or more, then G consists of 

the two straight segments xp and fry. 

(b) Otherwise, G is the preimage of the straight line from f{x) to 
f{y), where f is the quadratic map in the plane given by x,y,p 
to the Riemann Surface. 

5.2 Discretizing a Path 

First we approximate any (x, y)-path 7 with a piecewise linear (x, y)-path 
a that is composed of segments that are sufficiently long. The properties of 
a, formalized in the following definition, are used in the rest of this section 
to obtain a piecewise linear approximation using Steiner points. 

Definition 5.1. Let P = {pi,P 2 , ■ ■ ■ jPn} be a set of points in and let 
x,y G P. Assume j is a path from x to y that is internally disjoint from P. 
Also, let 0 < Ea < I be a real number. A sequence {0 = to) ^ tk, tfc+i = 
1 } is a e^-discretizing sequence 0/7 if it has the following properties: 
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1. The nearest neighbor of both 7 (ti) and 7 (^ 2 ) is x, while the nearest 
neighbor of both ^{tk-i) and 7 (tfc) is y. 

2. For all 1 <i < k, we have d( 7 (tj), 7 (ij+i)) = £„ • N( 7 (tj)). 

In this case, the piecewise linear path a = {x, 7 (^ 1 ),..., 7 (ifc), y} is a Sa- 
discretized path 0 / 7 . 

The following lemma guarantees the existence of ea-discretizing sequences 
and, hence, £Q-discretized paths. Its proof is essentially the same as the proof 
of Lemma 4.2. 


Lemma 5.2. Let P = {pi,P2,-'' :Pn} be a set of points in and let 
x,y G P. Let 7 6 e a path from x to y that is internally disjoint from P and 
0 < £« < 1 a real number. There exists a £ a-discretizing sequence of 'y. 


Next, we show that the N-length of a £ci-discretized path of 7 is not 
much longer than the N-length of 7 . 


Lemma 5.3. If a is a ea-discretized path of j for some Sa < 1/2, then 

^N(a) < (1 + 4£q,)%(7). 

Proof. For 1 < i < A; — 1, let Oj be the line segment ( 7 (L), 7 (^+ 1 )). Fur¬ 
thermore, let ao and be the line segments (x, 7 (ti)) and {'y{tk),y), re¬ 
spectively. Also, let 7 i = 7 [tj,L+i], for 1 < i < — 1 , and 70 = 7 [ 0 , ti] and 

7fc = 7[4, !]■ 

We prove the stronger statement that %(«,) < (1 -|- 4£Q,)I'N(7i), for any 
0 < i < k. For i G {0,A;}, the straight line is indeed the shortest path in 
the nearest neighbor metric (since N( 7 (ti)) = 7 (to) = x and N{'y(tk)) = 
'y{tk+i) = y), and so, I'N(«i) < ^N(7i)- 1 < ^ < A: — 1, by Lemma 5.4, 


^N(aj) < iNili) ■ 


N( 7 (L)) P£{ai) 
- (.{aiY 


Then, by the second property of £Q-discretizing sequences, 


^N(ai) < ■^N(7j) • 


N(7(A0) + £aN(7(L)) 
N(7(Ai)) -eaN(7(ti))' 


After simplifying and using the fact that £q, < 1/2, we obtain 


%(«*) < .^N(7j) • 


1 + £« 
1 £a 


< ■^N(7i) • (1 + 4e«). 


□ 
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5.3 Approximating with Steiner Points 

Assume P C x,y € P, and let 7 be an arbitrary (x,y)-path. We show 
how to approximate 7 with a piecewise linear path through a collection of 
Steiner points in To obtain an accurate estimation of 7 , we require the 
Steiner points to be sufficiently dense. The following definition formalizes 
this density with a parameter 5. 

Definition 5.2 (5-sample). Let P = {pi,P2r'' iPn} be a set of points in 
and let D C For a real number 0 < 5 < 1, a 5-sample is a (possibly 
infinite) set of points TOD sueh that if z G D\P, then d(z,T) < 5 -N( 2 ;). 

The following lemmas guarantees that an accurate estimation of 7 can 
be computed using a 5-sample. 

Lemma 5.4. Let P = {pi,P2r" ^Pn} be a set of points in and let 
x,y G Let 7 be any path from x to y, and I be the straight {x, y)-segment. 
Assume further that N(x) > i{l). Then, 


£n(0 < ■^n(7) • 


N(x)+£(0 
n(x5 - 


Proof. Let 7' = 7 n B{x,i{l)). Observe that 7' is a finite collection paths, 
and let ^( 7 ') be the total length of these paths. For any point z G B{x, l{l)) 



Figure 4: Proof of Lemma 5.4; 7 ' is bold. 

we have N(x) — (.{1) < N(z) < N(x) -|- i{l) because N(-) is Lipschitz. Since 
and I are both contained in B{x,i{l)), in particular, we have: 

%(7')>^(7')(N(x)-^(/)) 

and 

^n (0 < ^( 0 (N(x) + m) < %')(N(4:) + i{l)). 
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Thus, since 7 ' C 7 , we have •^n(70 ^ ^n( 7 )- Finally, because both sides of 
the inequalities are positive numbers, we have: 


^n( 0 ^ %(0 ^ N(x)+£(0 
%(7) “ ^n( 7') “ N(x)-£(/)■ 


□ 

Lemma 5.5. Let P = {pi,P 2 , ■ ■ ■ ,Pn} be a set of points in and let S 
be a 6-sample, and let 0 < 6 < 1/10. Then, for any pair of points x,y & P, 
there is a piecewise linear path 7 = (x, si,.. ., Sk,y), where si,. .., Sfc G S, 
such that: 

^n( 7) < (l + Ci52/3)dN(x,y), 
and, for all 1 < i < k — 1, 

£]<i{{si, Si+i)) < C 2 ■ 6"^/'^ ■ N(si). 


Cl and C 2 are universal constants. 

Proof. Let 7 be an (x, y)-shortest path under the nearest neighbor metric. 
We assume, without loss of generality, that 7 is internally disjoint from P. 
Of course, we can break any path to such pieces and use induction to prove 
the lemma for the general case. 

Let a = (x = 7(^0)) 7(^1)) • ■ ■ ilitk+i) = y) be a ( 2 J^/ 3 j_discretization of 
7. For all 0 < i < k, let a* = (7(^1),7(tj+i)) for 0 < i < A:, let m* be the 
midpoint of a*, and let Sj G S' be the closest Steiner point to m^. Then, let 
fii = {'y{ti), Sj, 7(tj+i)) be a two-segment path, and let /3 = / 3 o o / 3 i o • • • o jdj.. 
Finally, let 7 = (x, si, S2, ■ ■ ■, Sfc-i) v), and let pi = {si, Si+i) where 0 < i < k 
and So = X and = y to simplify notation. 

We need to prove that ^n( 7 ) < (1 + C'i 5 ^/^)£n( 7 ). Lemma 5.3 implies 
that 

^N(a) < (1 + 85 ^/^)£n(7) 

Showing the following two facts imply the statement of the lemma. 

^n(/?) < (l + 20<52/3)£N(a) (5) 

%(7) < (1 + 852/3)£n(/3) (6) 
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Figure 5: The lines segment Oi is blue, and the two-segment path Pi is red. 


A proof of Inequality (5): We prove the stronger statement that for all 
0 < t < /c, ^n(A) < (1 + Recall that m* is the midpoint of Oj 

and Sj G 5 is the closest Steiner point to m-j. Let hi be the closes point of 
ai to Si] see Figure 5. 

Because l{ai) k. N(mj), i{{mi, Si)) 6 N(mj), and 6 < 5^/^, we have 
(sj, hi) is orthogonal to Oj. Then, the following bound on the length of the 
piecewise linear path Pi = ( 7 (L), Sj, 7 (ti+i)) is implied by the Pythagorean 
theorem using derivative to maximize a function. 

m) < 2 • < 2 • Y^(<52/3N(7(t,)))' + (5N(m,))' (7) 

Then, the Lipschitz property of the N-function implies: 

NK) < N(7(ti)) + ^ < N(7(L)) + 52/3 . N(7(L)) < (1 + <52/3)N(7(ti)). 

By substituting (1 -|- 6‘^^^)N{'y{ti)) for N(mi) in (7), we obtain 

i{Pi) < 2 • (h2/3N(7(L)))" + (,5(1 + <52/3)N(7(ti)))' 

= <52/3n( 7(L)) • y^l + (hi/3(l + 52/3))2 

= £iai) • (1 + 252 / 3 ) 


On the other hand, because both Oj and Pi are contained in B{'^{ti),ea^{'p{ti))), 
for any point z G aj U /3j we have N(7(tj))(l — 2 ^ 2 / 3 ) < N(2:) < N(7(tj))(l -|- 
252/3) and so 

^n(A) <W-N(7(ii))(l+ 252 / 3 ). 

and 

%(ai) >^(ai)-N( 7 (ti))(l-252/3). 
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Last three inequalities together imply 

Wi) < l + / 

%(«*) - 1 - 252/3 ■ V 

< (1 + 205^/^). 

Thus, we obtain Inequality (5) by summing over all a^’s and /3j’s. 

A proof of Inequality (6): First, by the triangle inequality, we obtain: 

< iisi,ti+i) + e{ti+i,Si+i) 

N(7(L)) + N(7(L+i)) ^ 

< 25^/^ • max(N( 7 (ti)),N( 7 (L+i))) 

Then, Lemma 5.4, letting x be the endpoint with larger N value, implies 

1 + 252/3 

■^N (.'Hi) — '^N('Si, L+l, •Si+l) ■ 252/3 ~ /*+l! 'Si+l) ' (1 T 85 ^ ). 

Finally, by summing over all ry^’s we obtain Inequality 6, and the proof is 
complete. □ 

5.4 The Approximation Graph 

So far we have shown that any shortest path can be approximated using a 
5-sample that is composed of infinitely many points. In addition, we know 
how to compute the exact N-distance between any pair of points if they 
reside in the same Voronoi cell of Vor(P). Here, we combine these two ideas 
to be able to approximate any shortest path using only a finite number of 
Steiner points. The high-level idea is to use the Steiner point approximation 
while 7 passes through regions that are far from P and switch to the exact 
distance computation as soon as 7 is sufficiently close to one of the points 
in P. 

Let P = {pi,P 2 , ■ ■ ■ ,Pn} be a set of points in and let B be any convex 
body that contains P. Fix 5 G (0,1), and for any 1 < 1 < n, let r* = rp{pi) 
be the inradius of the Voronoi cell with site pi. Also, let Ui = {1 — 

Finally, let S' be a 5-sample on the domain B \ Ui<i<n B{pi, Ui). 

Definition 5.3 (Approximation Graph). The approximation graph A = 
A{P,{ui,...,Un},S,5) = (TajL-^i) is a weighted undireeted graph, with 
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weight function w : E_a —)■ M'*'. The vertices in are in one to one corre¬ 
spondence with the points in S U P; for simplicity we use the same notation 
to refer to corresponding elements in S'UP and V^. The set P _4 is composed 
of three types of edges: 

1 . If si,S 2 e S and si,S 2 G B{pi,ri) for any pi, then (si,S 2 ) G Pyi and 
'w{si,S 2 ) = dN(si,S 2 )- We compute this distance using Corollary 1. 

2. Otherwise, ifsi,S 2 G S and i{si,S 2 ) < (72(5^/^ max(N(si), N(s 2 )), 
where C 2 is the constant of Lemma 5.5, then (si, S 2 ) G Pyi and w(si, S 2 ) = 
max(N(si),N(s 2 )) • ■^(si,S 2 )- 

5. If Si G S and si G B{pi,ri) then {pi,si) G Eji and w{pi,si) = 
dN(Pi,si) = (d(pi,si))2/2; see Corollary 1. 

For x,y gVji let dji{x, y) denote the length of the shortest path from x to y 
in the graph A. 

The following lemma guarantees that the shortest paths in the approxi¬ 
mation graph are sufficiently accurate estimations. 

Lemma 5.6. Let {ui ,..., Un}, S and 5 he defined as above. Let A{P, {ni,..., Un}, S, 5) 
be the approximation graph for P. For any pair of points x,y G P we have: 

(1 - (72(5^/^) • dN(x, y) < d^(x, y) < {I + ■ dN(x, y), 

where C 2 and C 4 are constants computable in 0(1) time. 

Proof. First, we prove the upper bound. Let D = B \ P(pj, n*), 

and recall that S is a (5-sample on D. Let S+ be a 5-sample of B such 
that S"*" D D = S. According to Lemma 5.5, there is a piecewise linear 
path y = {x, si,..., Sk,y}, for si,...,Sfc G S~^, such that ^n(^) < (1 + 
(7i5^/^)dN(x, ?/). We use y to obtain a (x,y)-path a in A with length at 
most (1 -|- 2(725^/^)1 'n(^); O 2 is the constant in Lemma 5.5. 

Consider any maximal subpath 7 /' = (si, Si+i,..., Sj-i, Sj) of y that re¬ 
sides in B {pt, Ut/{1 — 5^/^)) for a pt G P. The lengths of the type (2) edges 
and the second inequality of Lemma 5.5 imply that Si, Sj ^ B{pt,ut), which 
in turn implies that Si, Sj G D. Further, {si,Sj) is a type (1) edge of A 
and w{si, Sj) = Sj) < i^iy'). So, y' in y corresponds to a single type 

(1) edge {si,Sj) in a with no longer length. In fact, all type (1) edges of a 
correspond to not shorter paths in y. 

Similarly, we consider maximal paths {p, Si, Sj+i,..., Sj_i, Sj), where p 
is an input point, and replacing them with type (3) edges of a. Again, all 
type (3) edges of a correspond to not shorter paths of y. 
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It remains to show that type (2) edges of a are not too long. Observe, 
that each type (2) edge of a corresponds to a line segment of r]. Let (si, S 2 ) G 
be a type ( 2 ) edge and let max(N(si),N(s 2 )) = N(si) without loss of 
generality. By the Lipschitz property and the length property of type (2) 
edges 

^n(si,S2) > (N(si)-^(si,S2))^(si,S 2) > (l-C'2<5^/^)N(si)^(si,S 2 ) = (l-C2(5^/^)rc(si,S2) 
which implies 

S 2 ) < (1 + 202 ( 5 ^/^)%(si,S 2 ), 

and so, by Lemma 5.5, 

iA{^) < (1 + 2 C 2 < 52 / 3 )£N(r?) < (1 + C' 4 < 52 / 3 )dN(x, y). 

To prove the lower bound, let x,y £ V_a_ and let a = {x, si,..., Sk,y) 
be a path from x to y in A. Also, let 7 be the (x, y) path in that hits 
(x = So, si,..., Sfc, y = Sfc+i) in this order, and let the {si,Sj) subpath of 
7 be a line segment if {si,Sj) is a type ( 2 ) edge in A and be the shortest 
{si,Sj) nearest neighbor path otherwise. Then, let 7 * be the subpath of 7 
from Si to Sj+i. We compare w{si, Sj+i) with I'N( 7 i). 

If (sj, Si+i) is type (1) or (3) then ^N( 7 i) = w{si, S 2 ). Otherwise, (sj, Sj+i) 
is type (2) and by Lipschitz property 

^n(si,S2) < (N(si)+£(si,S2))^(si,S2) < (1+C2(^^/^)N(si)£(si,S2) = {l+C25‘^/^)w{si,S2) 
Overall, 

dN(x, y) < £n(7) < (1 + C 2 d^^^)£A{<^), 
and the proof is complete. □ 

5.5 Construction of Steiner Points 

The only remaining piece that we need to obtain an approximation scheme 
is an algorithm for computing a (5-sample. For this section, given a point set 
T and x G T, let rr(x) denote the inradius of the Voronoi cell of Vor(r) that 
contains x. Also, given a set T and an arbitrary point x (not necessarily in 
T), let ^t{x) denote the distance from x to its second nearest neighbor in 
T. 

We can apply existing algorithms for generating meshes and well-spaced 
points to compute a (5-sample on V\\^- B{pi, Ui), where P C is a domain, 
and Ui = {1 — 6‘^^^)rp{pi). The procedure consists of two steps: 
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1. Use the algorithm of [MSV13] to construct a well-spaced point set M 
(along with its associated approximate Delaunay graph) with aspect 
ratio r in time (nlogn-|- \M\). 

2. Then over-refine M to S for the sizing function g{x) = j^fp(x) (while 
maintaining aspect ratio r) in time 2 ^*^'^)|S'| by using the algorithm of 
Section 3.7 in [Shell], (see also [HOMS 10] for an earlier use of this 
technique) 


In the above algorithm, we will choose r to be a hxed constant, say, r = 6 . 
Both of the meshing algorithms listed above are chosen for their theoretical 
guarantees on running time. In practice, one could use any quality Delau¬ 
nay meshing algorithm, popular choices include Triangle [She96] in and 
Tetgen [Sill] or COAL [ART+12] in 

From the guarantees in ([Shell]), we know that 

where A is the spread of P, i.e., the ratio of the largest distance between 
two points in P to the smallest distance between two points in P. 

Lemma 5.7. If x £V \ B{pi, Ui), where pi = NN(x), then N(x) < fp{x) < 
5N(x). 


Proof. Note that N{x) < fp{x) is trivial. To prove the second half of the 
inequality, let q be the closest point to j? in P \ {pi}. Then, by the triangle 
inequality, we have that 


fp{x) 


< 

d(x, 

<q) 


< 

d{x. 

,Pi) + 

d(K,g) 

< 

d{x, 

<Pi) + 

‘2rp{pi) 

< 

d{x, 

<Pi) + 

2 

1 - 52/3 

= 

5d{x,Pi) 


= 

5N(: 

x). 



□ 

Now, it remains to show that the point set S is indeed a 5-sample on 
V\yj^B{pi,Ui). This is provided by the following lemma. 

Lemma 5.8. S is a 5-sample on D \ jj■ B{pi, m). 
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Proof. Let x G V \ \Ji B{pi,Ui) be an arbitrary point. Then, let q be the 
closest point to x that lies in S. Note that 


d{x,q) < T-rs{q) 


< T 





This implies that 


d{x,q) < 


< 



< 


(5N(x), 


as desired. 


□ 


Now, we calculate the number of edges that will be present in the ap¬ 
proximation graph defined in the previous section. For this, we require a 
few lemmas. 

Lemma 5.9. Let A = B{pi,rp{pi)) \ B{pi,Ui) be an annulus around pi. 
Then, |An5| = 

Proof. By the meshing guarantees of [Shell], we know that for any point 
s G A n S', B{s, t) does not contain a point from S \ {s} for t = Q{rs{s)) = 
• rp{p)). Thus, the desired result follows using a simple sphere packing 
argument. □ 

Lemma 5.10. If s G S, then \B{s,C 26 ^/^N{s)) n S| = where C 2 is 

the constant in Lemma 5.5. 

Proof. As in the previous lemma, meshing guarantees tell us that for any 
s' G B{s,C 2 d^/^N{s)), we have that B{s',t) does not contain a point from 
S\ {s'} for t = n((5 • N(s')) = n((5 • N(s)). Thus, we again obtain the desired 
result from a sphere packing argument. □ 

From the above lemmas, we see that A is composed of jSj = A 

vertices and + |5'| • = |5'| • edges. 

Remark. Note that the right hand side of (8) is in terms of the spread, 
a non-combinatorial quantity. Indeed, one can construct examples of P for 
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which the integral in (8) is not bounded from above by any function of n. 
However, for many classes of inputs, one can obtain a tighter analysis. In 
particular, if P satisfies a property known as well-paced, one can show that 
the resulting set S will satisfy |5| = (see [MPS08, Shel2]). 

In a more general setting (without requiring that P is well-paced), one 
can modify the algorithms to produce output in the form of a hierarchical 
mesh [MPSll]. This then produces an output of size and (1 -|- e)- 

approximation algorithm for the nearest neighbor metric can be suitably 
modified so that the underlying approximation graph uses a hierarchical set 
of points instead of a full d-sample. However, we ignore the details here for 
the sake of simplicity of exposition. 

The above remark, along with the edge count of A and the running time 
guarantees from [MSVI3], yields Theorem 1.2, the main theorem of this 
section. 


6 Discussion 

Motivated by estimating geodesic distances within subsets of M"', we consider 
two distance metrics in this paper: the V-distance and the edge-squared 
distance. The main focus of this paper is to find an approximation of the 
V-distance. One possible drawback of our (1 -|- e)-approximation algorithm 
is its exponential dependency on d. To alleviate this dependency a natural 
approach is using a Johnson-Lindenstrauss type projection. Thereby, we 
would like to ask which properties are preserved under random projections 
such as those in Johnson-Lindenstrauss transforms. 

We are currently working on implementing the approximation algorithm 
presented in Section 4. We hope to show that this approximation is fast in 
practice as well as in theory. 
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